Pressure dependence of the melting mechanism at the limit of overheating in 

Lennard-Jones crystals 
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We study the pressure dependence of the melting mechanism of a surface free Lennard-Jones 
crystal by constant pressure Monte Carlo simulation. The difference between the overheating 
temperature(Tbff) and the thermodynamical melting point(TAf) increase for increasing pressure. 
When particles move into the repulsive part of the potential the properties at Toh change. There is 
a crossover pressure where the volume jump becomes pressure-independent. The overheating limit 
is pre-announced by thermal excitation of big clusters of defects. The temperature zone where the 
system is dominated by these big clusters of defects increases with increasing pressure. Beyond 
the crossover pressure we find that excitation of defects and clusters of them start at the same 
temperature scale related with Toh- 

PACS numbers: 64.70.Dv, 61.72.Ji, 65.40.-b 
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I. INTRODUCTION 

Melting is a dramatic phase transition where the trans- 
lation symmetry present in the solid phase is destroyed at 
the transition point. As a paradigm of a phase transition 
in condensed matter physics, its microscopic mechanism 
has been studied from the beginning of the solid state 
theorjii. 

Early theories relate melting with an intrinsic elastic 
instability of the crystal. According to a stability crite- 
ria established by Bom 2 , melting takes place when an 
elastic shear modulus vanishes and the crystal looses its 
ability to resist shear. This was interpreted as due to the 
softening of a transverse phonon mode with, as a con- 
sequence, the homogeneous breakdown of the crystalline 
order. This notion of crystalline instability refers only 
to the properties of the crystalline phase but a complete 
theory of melting should also refer to the properties of 
the liquid phase. 

Therefore, modern theories invoke thermal excitation 
of extended defects such as dislocation lines^. Pairs 
of dislocations can be thermally excited at temperatures 
within the range of typical solid phases. The presence 
of dislocations lowers the cost in energy to create an 
additional dislocation; if the energy reduction is strong 
enough, it could lead to an avalanche of dislocations re- 
sulting in a first order phase transition. 

So far, two different ways of melting has been recog- 
nized in the literature^. One related to the thermody- 
namical melting starting at the surface or extrinsic in- 
homogeneities of the samples, and the other one with 
regards to the limit of the mechanical stability of the 
crystal which takes place when surface is avoided and the 
solid could be overheated above the equilibrium melting 
point. Overheating has been experimentally achieved by 
coating a solid with another material of higher melting 
point 6 . 

As no extrinsic inhomogeneities are present, thermal 
excitations of defects play an essential role to initiate 



the melting process. In this work, we will focus on this 
type of melting. 

Although experimental analysis of the microscopic dy- 
namics near melting is very difficult, numerical simula- 
tions have recently shown that correlated clusters of de- 
fects thermally excited play a central role in this process 
at the limit of overheating-. Moreover in a previous work 
we have related clusters of defects to dislocation lines as- 
sumed in phenomenological theories^. However, the sit- 
uation is far from being solved and continues to puzzle 
condensed-matter theorists. 

An alternative strategy to look for the melting origin 
could be to change an external parameter and to fol- 
low how this parameter affects the melting. The obvious 
choice of this parameter is the external pressure acting 
on the solid. Moreover, the pressure dependence of the 
melting temperature has a practical interest per se and 
the determination of high pressure melting curves is a 
subject of many investigation al' 1 ? . 

The purpose of this paper is to discuss the effect of 
pressure on the melting properties of a surface free crystal 
by means of a computationally intensive Monte Carlo 
(MC) calculation using a constant pressure algorithm. 

The paper is organized as follows. In Sec. II we analyze 
some thermodynamic properties at overheating limit. In 
Sec. Ill we look at the defect structure preceding the 
overheating limit. Section IV contains our conclusions 
and discussion. 



II. PHASE DIAGRAM OF THE 
LENNARD-JONES SYSTEM 

Our simulations have been performed on a cubic box 
of 2048 particles 11 interacting via a Lennard-Jones (LJ) 
potential written as V(r) = 4e[(2) 12 - (f ) 6 ]. To guar- 
antee that the potential goes smoothly to zero at dis- 
tance greater than the cutoff (r c ), a cutoff region from 
0.95r c to r c was used following previous work<A2*A£. There- 



fore, we can compare our results for the properties at the 
overheating limit with ones previously obtained at the 
equilibrium melting point. Moreover we fix the cutoff 
r c = 2.1cr which has been shown to be greater enough to 
assure convergence of the results^. 

The energies are measured in units of e, the distances 
in units of V4cr and the pressure in unit of e/4cr 3 . The 
temperature is reported in units of e/feg. As a con- 
crete example we can take the parameters for argon 
(a = 3.4 A, e = 0.0104eV) where our units are: 5.397A 
for the distances, 120. 64K for the temperatures and 
10.5 x 10 6 Pa for the pressures. Our units of pressure dif- 
fer by a factor j from those used in previous worksi2ii^. 
Moreover, our unit of volume differ by a factor 4. 

We emphasized that the use of periodic boundary con- 
ditions (PBC) suppresses surface effects and the phase 
transition temperature we find does not correspond to 
the thermodynamical melting point but to the limit of 
overheating of the crystal^. 

Let us discuss how we constructed the phase diagram of 
the system. For a fixed pressure at a given temperature, 
we have averaged the internal energy and the volume over 
the MC runs. Both the particle coordinates and the vol- 
ume of the cell are updated in each MC step^. Heating 
is done in a step-by-step procedure starting at low tem- 
perature with the particles at the sites of a fee perfect 
lattice. The calculation of the jump for different proper- 
ties (volume, energy, etc.) at the limit of overheating is a 
difficult task in MC simulation, because of the long CPU 
time necessary to achieve equilibrium values in this in- 
termediate region. We have densely partitioned the tem- 
perature interval and run ~ 10 5 MC steps in this zone. 
In Fig. da) we show a typical results of the evolution 
of the energy and the atomic volume (v) for P = 1000 
with the temperature. The abrupt jump is related to 
the limit of overheating. We have observed that both 
the snapshot of the particle positions as well as the ra- 
dial distribution function (RDF) above this jump do not 
show any indication of the crystalline order. By locating 
the temperature where this jump takes place at different 
pressures we construct the figure ^b) where we show 
the overheating limit {Ton) as a function of pressure. For 
comparison, we have included the results for the equilib- 
rium melting point (Tm) extracted from Table II of Rcf. 
IT3I obtained by using molecular dynamics (MD) coexis- 
tence simulations. We take from this work the results 
obtained using the same smooth cutoff r c = 2.1a as in 
our MC simulation. 

The two set of values could be fit following the corre- 
sponding next laws: 

T M = 0.036 P 805 + 0.66 (1) 

Toh = 0.044 P 806 + 0.78 (2) 

Note in figure D^b) that even at P = 0, the overheating 
phenomena is observed. While the difference between 
Tm and Tqh increase with pressure, the exponent of 
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FIG. 1: (a) atomic volume (square, left scale) and internal 
energy (circles, right scale) vs. the temperature (T) for pres- 
sure P = 1000 (b) the overheating temperature (Toh) as a 
function of the pressure . We also include the melting point 
Tm as obtained in Ref. ^3 (open diamond) and from Ref. ^3 
(starts). Dashed (dot-dashed) lines correspond to the power 
law fitting given in eq. (1) (eq. (2)) shown in the text. 

the pressure is quite similar, ~ |. In the same fi gure 
we have added Tm datas (start symbols) from Ref. Ha . 
where the Gibbs-Duhem integration molecular simula- 
tion thechnique was used. We can see the agreement 
with the extrapolation of eq. QJ even in the high pres- 
sure regime^. 

Let us analyze this exponent at the light of previous 
results valid in the regime of high pressures*. When par- 
ticle are close enough, the repulsive part of the poten- 
tial ^2 dominates and a power law has been deduced 
(P/?4 = Const). All this, indicate that the low pressure 
zone could be neglected in the scale we are presenting the 
phase diagram. 

A microscopic understanding of the differents multi- 
plicative constant shown in the eq. and eq. |2"| 
could be explained with the following argument. While 
no extrinsic nucleation center as surface or defects are 
present in our simulation system, in order to break the 
crystalline order thermal activation of clusters of defects 
are necessary (see next Section). However, as pressure 
increase, more thermal energy is needed because parti- 
cles are strongly compacted. This is different than the 
process involved in the thermodynamical melting where 
a liquid phase starts at the surface and pressure is less 
effective to increase the transition temperature. 

The Born criteria has been recently reobtained for ho- 
mogeneous lattices under an arbitrary but uniform ex- 
ternal loacUii^. In the case of hydrostatic pressure the 
generalized stability criteria are given by: 

B T = (Cii + 2C7i 2 + P)/3 > (3) 
G = C 44 - P > (4) 
G' = (C11 - Ci 2 - 2P)/2 > (5) 

where Cn, C12 and C44 are the three different elastic 
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FIG. 2: (a) The atomic volume of the solid (voh) at the 
overheating limit as a function of P (dashed line, square sym- 
bols). The solid line shows the volume where G' vanishes at 
T = 0. We also include the volume at the melting point (vm) 
extracted from Ref. |l2 (dot dashed line, diamonds). The in- 
set shows the volume dependence of the cohesive energy of an 
ideal fee crystal. The dashed (solid) arrow signal the thermal 
energy to melt (overheat) the crystal, (b) The percentage 
jump of the atomic volume at the overheating limit. 



constants. The Bt, G and G are the generalized bulk 
modulus and the two different shear modulus. We have 
obtained the elastic constants by means of the harmonic 
approximation — on the Lennard- Jones fee crystal with 
interactions up to the third nearest neighbour and vari- 
able v. This procedure gives the elastic constants at 
T = under pure dilation or compression. Note that 
for a general value of v, the crystal is not in equilibrium 
under the action of the interacting forces and the Cauchy 
relation Cn — C44 is not in general satisfied. For each 
pressure we look for the value off where some of the mod- 
ule of Eq.(|3l5p vanish. At P = 0, B T is the first to be 
vanished, i.e. for smaller v, than the two shear modules. 
This instability take place at v = 0.315 higher than the 
melting detected by MC simulation (v ~ 0.273). From 
P > 15 the shear modulus G is the first to become zero. 
This v signals the volume at T — where the system 
could not resist shears. In Fig|2fa) we show these val- 
ues (solid line) together with the value of v correspond- 
ing to the overheating (voh) obtained by MC simulation 
(square symbols and dashed line). 

We have also included in the figure the volume at the 
melting point (t;m)- The density ( p)-temperature (Tm ) re- 
lationship were extracted from Ref. [l2j. The specific 
volume is v — - and the pressures at melt were obtained 
from Tm using eq. (Q. 

The three curves intersect at a pressure which sepa- 
rates two regimes. At this pressure, P c ~ 100, voh is 
approximately the value where a perfect fee lattice is in 
equilibrium under the action of the interacting forces. 
Therefore the system is dilated at melting for P < P c . 
For P > P c , it is compressed. Note that the stability cri- 
teria given by Eq.JSJ) gives an overall correct estimation 



for the overheating volume. The fact that the vanishing 
of G shear module gives a good estimation for the super- 
heating volume when solid expands has been recognized 
in previous worksLi&SI. Note however that in a recent 
study of amorphization under decompression it has been 
shown that the condition G' = and the softening of a 
shear phonon signal the critical pressure which destroys 
the crystalline ordesSi, only at T = 0. 

Note in addition, that vm is smaller (greater) than voh 
for the dilated (compressed) system. This is consistent 
with the overheating phenomena found for all pressure. 
Looking at the inset it is easy to note two different be- 
haviors around the equilibrium volume (v e ). For volumes 
bellow the v e the thermal energies necessary to destroy 
a crystal are lower in the overheating case. And in con- 
trary, for volumes above v e the thermal energies are lower 
in the melting case (see the inset of fig. 

The crossover from dilation to compression is also seen 
in the percentage jump of the atomic volume at the melt- 
ing Av/voh{Av is the difference between the liquid and 
solid volumes per particle). In Fig. Efb) we show this 
value obtained from MC simulation. We can see that for 
P>>P C , Av/voh is reduced and becomes almost pres- 
sure independent. This behavior could be understood 
from the relation valid in the high pressure regime. From 
Eq. (18) of Ref. [l3| we see that the specific volume of 
the solid (vs) and the liquid (vl) scale with the same 
exponent of the temperature: 

Av L = dpi 

4v s = C 2 pi (6) 

with $ = C x = 1.23 and C 2 = 1.18. The relative 
jump of the volume is: 



which is T and P independent and of the order of mag- 
nitude of the value found in the MC simulation at high 
pressure as can be seen in Fig. |5] (b) 

Microscopically, the structure of the solid near the 
overheating limit is not the same above and below P c . 
Indeed, the jump percentage of v becomes very small 
above P c signalling that the structure of the solid below 
Toh has a lot of defects and could be easily destroyed. 
We will return to this point in the next section where the 
evolution of the defect structure will be analyzed. 

We now look at the two term of the Claussius- 
Clapeyron equation obtained from our numerical results. 
Note that the validity of this equation is not guaran- 
teed a priori because we are calculating the properties at 
the overheating limit where the free energy of the liquid 
and solid phase are not equal. It is an interesting point 
to analyze how different arc the properties at Toh from 
the ones predicted assuming coexistence between the two 
phases. We have calculated the latent heat of the transi- 
tion defined as Lh = Aw + PAw, where Au is the jump of 
the internal energy. In Fig. [3] we compare the two terms 
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FIG. 3: The two terms of the Clausius-Clapeyron equation 
evaluated numerically. The solid line corresponds to 
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evaluated from Eq. [5] Squares correspond to T °"^ v where 
the latent heat Lh is obtained from the jumps of the internal 
energy and the volume (Aw). 



with ^m. 

dP 



of the Clausius-Clapeyron equation TAv/L^ 
obtained by derivation of the fitting law Eq. Striking 
the two term coincide at high pressure. 



III. STRUCTURE OF DEFECTS PRECEDING 
THE OVERHEATING LIMIT 

We have recently found^ a premelting temperature at 
P = where the number of thermally activated defects 
increases dramatically and defects start to group in clus- 
ters. Near Tqh only one cluster across the whole system 
survives, breaking the crystalline order. In this work we 
analyze how this feature evolves with the external pres- 
sure. 

We define a defect as a particle with coordination num- 
ber (CN) different from 12, which is the number of near- 
est neighbors (NN) in an ideal fee lattice. This CN, called 
Cnn, is obtained by counting the number of particles 
around a given particle up to a cutoff radius. This cutoff 
is chosen for each pressure as the value where the ra- 
dial distribution function has its first minimum at low 
temperature, i.e. the value of the distance between the 
maxima of the first and the second neighbors. 

To compare our results for different pressures we nor- 
malize the temperature scale to the value of Ton corre- 
sponding to each pressure (T = T/Tqh)- In Fig. 0] we 
show the coordination number mean value and the per- 
centage of defects as a function of T for different pres- 
sures. We see that the reduced temperature where a 
great number of defects are activated and the crystalline 
order is perturbed, decreases by increasing pressure. For 
example for P = we should go up to 70% of Ton to 
have 10% of defects whereas at P > P c we achieved the 
same quantity of defects by going only up to 50% of Ton- 
The decrease of the coordination number shows a simi- 
lar behavior. These facts are consistent with the small 




FIG. 4: The percentage of defects (left scale) and the mean 
coordination number (right scale) as functions of the reduced 
temperature for different pressures. 



jump of the volume between the two phases found in the 
previous section. 

Moreover, it is seen in Fig. 0] that going beyond the 
crossover pressure P c the curves almost coincide. When 
the repulsive part of the potential dominates the system 
could be connected with a hard sphere model with a tem- 
perature dependent radius22. As we normalize the tem- 
perature scale at the Tor, the coincidence of the defect 
structure at high pressure is a manifestation of the equiv- 
alence of the models to an effective hard sphere model. 

Like we stated in the introduction, modern theories 
associate the melting process with a proliferation of dis- 
location lines. This means that defects would have to 
be correlated near the melting point. To analyze this 
correlation, we have grouped the defects into clusters by 
the following methodology. We start from a given de- 
fect and search for new defects up to the cutoff distance. 
For each of these new defects, we undertake the same 
procedure. We iterate this process up to completion of 
a cluster of connected defects. Then we take a new de- 
fect disconnected to all the previous ones and develop the 
same procedure. At the end of this process we separate 
our set of defects in N c i clusters which arc disconnected 
between them. The defects within a cluster are neighbors 
to each other but they are not connected to the defects 
of another cluster. 

In Fig. [3 we show the mean value of N c i as a function of 
T. The decrease of this quantity above a given T (called 
T pm ) should be compared with the increase of the num- 
ber of defective particles shown in Fig. These facts 
indicate that the clusters are becoming bigger and that 
the defects correlate among them for T > T pm . T pm de- 
creases with increasing pressure and becomes constant at 
high enough pressures. The premelting zone increases as 
we already emphasized previously. This could be impor- 
tant to detect experimentally excitations of dislocations 
lines in the study of melting of materials under pressure. 

For low pressures we see that the curves move away 
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FIG. 5: The average number of clusters as a function of T 
for different pressures. 

from the P = curves and for high pressures they no- 
tably coincide with each other. All of these seem to point 
out that the same mechanism makes the crystal melt at 
these values of pressure. Note that in a phenomenological 
theory the crystal melts at a fixed density of dislocation 
lines at different pressures^. At this stage, we cannot 
give a definitive statement on this prediction but it is 
suggestive that the tails of all the curves of Fig. ^corre- 
spond to only one cluster. This means that all the defects 
of the system are interconnected to each other. 

Finally, let us notice that Fig. 0] and Fig. [S] show not 
only the same mechanism of melting at different pressures 
but also the same structure of thermal excited defects at 
these pressures. 

IV. CONCLUSIONS AND DISCUSION 

We have studied the effect of pressure on the overheat- 
ing properties of a crystalline solid. We have compared 
the properties at the overheating limit with the ones at 
the thermodynamical melting where the solid and liquid 
could coexist. We have connected the thermodynamic 



properties with the microscopic dynamics near the over- 
heating limit. We have taken as a representative example 
an fee crystal whose particles interact by a Lennard- Jones 
potential. We have shown that it is possible to obtain 
equilibrium properties of both phases in the overheating 
zone by constant pressure MC simulation. 

Important conclusions of our study are: 

(a) The difference between Toh and Tm increase with 
increasing pressure. 

(b) A similar power law fit both the pressure dependence 
of Toh and T M , 

(c) Volume expansion due to thermal effects and volume 
compression due to the external pressure balance at a 
given pressure. This pressure signal a crossover between 
two regimes. The volume jump between the two phases 
becomes rather small(~ 5%)well above the crossover 
pressure. 

(d) The previous result is interpreted at microscopic 
level as due to the existence of some liquid-like struc- 
tures in the solid phase . These structures appear as 
defects in the crystalline phase. Indeed our results show 
that a great number of defects are activated above a 
prcmclting temperature. On the reduced temperature 
scale, this temperature decreases by increasing pressure 
and saturates at high pressures. 

(d) The prcmclting temperature also signals the be- 
ginning of a temperature range where defects correlate 
between them grouping in clusters. The number of 
clusters decreases as the overheating limit approaches 
from low T by cluster collapsing, leaving only one cluster 
near the overheating limit. 

Our results indicate that experimental detection of pre- 
melting behavior, which has been elusive up to now, 
could be easier detected in high-pressure studies of the 
melting. 
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